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ABSTRACT 

We present a comprehensive study of interstellar X-ray extinction using the extensive Chandra 
supernova remnant archive and use our results to refine the empirical relation between the hydrogen 
column density and optical extinction. In our analysis, we make use of the large, uniform data sample 
to assess various systematic uncertainties in the measurement of the interstellar X-ray absorption. 
Specifically, we address systematic uncertainties that originate from (i) the emission models used to 
fit supernova remnant spectra, (ii) the spatial variations within individual remnants, (Hi) the physical 
conditions of the remnant such as composition, temperature, and non-equilibrium regions, and (iv) 
the model used for the absorption of X-rays in the interstellar medium. Using a Bayesian framework 
to quantify these systematic uncertainties, and combining the resulting hydrogen column density 
measurements with the measurements of optical extinction toward the same remnants, we find the 
empirical relation A^h = (2.87±0.12) x 10^^ Ay cm“^, which is significantly higher than the previous 
measurements. 

Subject headings: ISM: dust, extinction - ISM: supernova remnants - X-rays: ISM 


1. INTRODUCTION 

The linear relationship between optical extinction 
(Ay) and hydrogen column density (A^h) has long been 
observed and utilized to estimate the X-ray or optical 
brightness for new sources or to fit the broadband spec¬ 
trum of X-ray sources. It is also used to obtain dis¬ 
tance estimates for X-ray sources using their measured 
column dens i ties ( s ee e.g..lDurant fc van Kerkwii3 |2006|: 
Giiver et al.l 20lA IRatti et al.l 120101 : ISoria et al.l 2012 : 
Nielsen et al. 1201^ . 

The photoelectric absorption by interstellar material 
causes rapid attenuation of soft X-rays in the spectra of 
Galactic sources. Measuring the extent of this attenua¬ 
tion yields information about the total column density 
along the line of sight to the source. Although this is 
commonly expressed in terms of the equivalent hydro¬ 
gen column density A^h, in the soft X-ray band (0.1 — 10 
keV), it is predominantly caused by abundant heavier 
elements such as O, Ne, Si, Mg and Fe. Optical extinc¬ 
tion is caused by grains of the same elements. Because 
dust tends to follow the metal distribution in the ISM, 
when averaging over many different lines of sight and 
over distances that are larger than the clumping scale 
of the ISM, it is reasonable to assume an approximate 
linear relationship between between Ay and A^h- 

There have been a number of different studies that 
have sought to accurately determine the relation between 
optical extinction and the hydrogen column density em¬ 
pirically. The methods vary across these studies, and 
the results show disc repancies greater than t he statis¬ 
tical errors for each. iReina fc Tarenghil (|1973ll used X- 
ray binaries and extended sources and found a relation 
of A^h = 1-85 X 10^^ X Ay cm~^ (hereafter, A^h is in 
units of cm“^ and Ay is in magnitudes); iGorenst^ 
(|1975ll used supernova remnants (SNRs) to find A^h = 


(2.22 ± 0.14) X 10^^ X Ay] while iPredehl fc Schmittl 
(1199511 used a combination of ROSAT point sources and 
SNRs and measured A^h = (1.79 ± 0.03) x 10^^ x Ay. 
Recently, iCiiver fc Ozell (12009( 1 collected a sample of 22 
SNRs, for which the hydrogen column density and op¬ 
tical extinction were previously measured, and found 
A^h = (2.21 ±0.09) X 10^^ X Ay. This had the advantage 
of using the high quality data from modern X-ray tele¬ 
scopes such as Chandra and XMM-Newton and focusing 
on sources with little-to-no intrinsic absorptiorQ, but had 
the disadvantage of being unable to account for system¬ 
atic errors that may vary for each published value. For 
example, when using values from the literature, there is 
no way to account for the variety of choices that are made 
during the data processing pipeline. The selected regions 
or the particular emission models may be interesting for 
the objectives of a particular study, but less ideal for 
the determination of the hydrogen column density and 
the comparison to the optical extinction measurements 
to determine the slope of the relation. 

In this study, we take advantage of the wealth of SNR 
data available in the Chandra archive to investigate and 
quantify the systematic errors present in the determina¬ 
tion of the hydrogen column density from the analysis 
of the X-ray data. We analyze all of the observations 
from the Chandra archive using standardized procedures 
so that we can quantify the systematic errors on each 
measurement. Using only observations from the Chandra 
archive also ensures a completely uniform data set: All 
observations were performed using the ACIS detector, 
and each had spectra generated and analyzed using the 

^ To illustrate this, consider a remnant with a 15 pc radius that 
has swept up 100 Msmi of ISM material (which is a large esti¬ 
mate for most SNRs). The column density through the remnant 
is approximately 2.6 X lO*^® cm~^, which is at least a factor of ten 
smaller than the A^h values presented in this paper. 
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Table 1 

Recovered A^h values from simulated data 


Assumed TVh 
( 10^2 cm-2) 

Ts 

(keV) 

Tes 

(keV) 

"^0 

(loll) 

Fit A^h 

(1022 cm-2) 

0.5 

0.3 

0.15 

10 

0.456 ± 0.094 




5 

0.412 ± 0.135 




1 

0.469 ± 0.052 


0.6 

0.3 

10 

0.525 ± 0.030 




5 

0.458 ± 0.023 




1 

0.499 ± 0.039 


0.9 

0.45 

10 

0.514 ± 0.020 




5 

0.492 ± 0.022 




1 

0.436 ± 0.032 

1.0 

0.3 

0.15 

10 

0.969 ± 0.139 




5 

0.922 ± 0.132 




1 

0.998 ± 0.143 


0.6 

0.3 

10 

1.007 ± 0.056 




5 

0.960 ± 0.056 




1 

0.988 ± 0.047 


0.9 

0.45 

10 

0.952 ± 0.051 




5 

0.973 ± 0.048 




1 

0.942 ± 0.061 


same treatment for background subtraction and model 
fitting routines using the spec tral analysis software xspec 
(version 12.8.1: lArnaudlfTQQGL with NEIVERS 1.1). This 
consistent treatment of a uniform data set gives us the 
opportunity to quantify the existing systematic errors in 
a way that had not been possible before. 

The uniformity of our data set also allows us to ex¬ 
plore the uncertainties associated with fitting spectra 
with a variety of models. SNRs can be home to a wide 
range of plasma properties due to the wide range of 
ages, host environments, and composition of the sam¬ 
pled gas. In general, however, SNRs exhibit continuum 
emission driven by thermal bremsstrahlung accompanied 
by emission lines produced by ejecta material or swept up 
ISM. In the 0.5-5.0 keV range, where we predominantly 
perform the spectral fits, these features manifest most 
visibly in magnesium, silicon and sulphur lines (though 
features due to iron, oxygen and neon are also possible 
within this range). Non-thermal features can also be 
present due to particle acceleration in SNR shocks. Due 
to the wide range in spectral features and plasma condi¬ 
tions, there also exist a range of models that can be used 
to fit X-ray spectra which will be detailed in Section 2. 

Making use of the high spatial resolution and large 
number of counts in the Chandra archival data, we also 
investigate the systematic uncertainties associated with 
differing lines of sight towards larger remnants. We as¬ 
sume that all absorption is due to interstellar gas and 
dust along a line of sight, with no intrinsic absorption 
from the optically thin remnant. This assumption allows 
there to be differing values for IVh for different regions 
in a given SNR due to real differences in gas along a 
line of sight, which is supporte d by deep observa tions 
of large remnants such as W49B (iSasaki e t al.ll20 I3ll and 
CTB 109 (also known as SNR G109.1-1.0: lKeohane et aH 
l2007f) . for which many regions can be fit. When the 
archival observations have sufficient duration to have 
multiple high-count regions, we are able to fit multiple 
regions to understand how the value of A^h varies over 
the remnant. By fitting multiple regions with a range 
of models, we can determine the magnitude of the total 



Figure 1. A Chandra image of SNR G109.T1.0 taken with the 
ACIS-I detector (data from OBSID 1901, red: 0.5-1.1 keV; green: 
1.1-1.6 keV; blue: 1.6-5.0 keV). The numbered circles show six 
of the 11 regions that were selected for spectrum extraction from 
this remnant. The dashed rectangle and circle shows the regions 
selected for background subtraction. 

systematic error on iVn measurements in the direction of 
the SNR. We note that any dust that is intrinsic to the 
SNR contributes negligible extinction. The integrated 
column density of swept-up dust is small, and that of 
dust condensed from the stellar ejecta is smaller still. 
Mor eover, SNR shocks are very efficient dust destroy¬ 
ers (jTemim et al.l[2M^ . making any contribution to the 
overall extinction exceedingly small. 

Einally, and possibly most importantly, many of the 
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Figure 2. Thermal spectra extracted from two different regions 
of SNR G109.1-1.0 (corresponding to regions 3 and 9 shown in 
Figure 0 in the 0.5 — 5.0 keV range with best fit model compo¬ 
nents overlaid. Region 3 spectrum exhibits weak emission lines 
(from the thermal component shown in the dotted blue line) with 
some power-law component contribution (shown in the red dashed 
line), while region 9 spectrum has stronger lines and no power-law 
component (thermal model shown as solid black line). 

previous analyses of C handra and other SNR data 
used solar abundances of I Anders fc Grevessd (|1989D de- 
rived from observ ations of the sun and meteorites. 
IWilms et al.l (j200r)fl showed that the abundances in the 
ISM can differ from solar values, which can affect the 
fit value of Nn. It is also important to consistently use 
an ISM absorption model that features improved calcu¬ 
lations of the absorption cross sections when seeking to 
improve systematic errors in the hydrogen column den¬ 
sity determination. We use the wilms abundance ta¬ 
ble and tbabs absorption model within xspec in order 
to accurately model the interstellar absorption to mea¬ 
sure the hydrogen column density. The tbabs model in 
xspec adopts the cross-sections from IWrner fc Yakovle'^l 
(jl995fl . This standard that we adopt for the entire sam¬ 
ple leads to larger 30%) column densities than those 
reported in the earlier studies. 

In Section [2l we discuss the variety of models used to 
fit our data, as well as the results of simulated data sets. 
Sections [3] and |4] detail the processes we used to process 
the data, as well as our approach to identifying system¬ 
atic errors and determining total uncertainties. Finally, 
in Section [S] we use our new data set to derive the slope 
of the Vh-Av relation. 

2. MODELS AND SIMULATED RESULTS 

The models we used to fit each thermal region 
{raymond, mekal, nei, and pshock) were chosen to 
span the range of complexity one might expect from 
a region of plasma in an SNR. The raymond model 
(|Ravmond fc Smitbl I1977D is the simplest, modeling a 
hot, diffuse gas in collisional equilibrium. Similarly, the 


mekal model also calculates spectra of a plasma in col¬ 
lisional equilib rium, but with imp roved handling of the 
Fe-L co mplex (iKaastra et al.lll996D . The nei and pshock 
models (jBorkowski et al.l l20nih are more complicated, 
providing spectra for a non-equilibrium ionization (NEI) 
collisional plasma (nei), and constant temperature NEI 
plasma heated by a plane-parallel shock (pshock). For 
non-thermal regions, the models powerlaw and/or srcut 
were used. Both models produce a spectrum we would 
expect from synchrotron emission from a power-law dis¬ 
tribution of shocked electrons interacting with an SNR’s 
magnetic field, with an additional expo nential cutoff of 
the distribution for the s rcut model (|Revnoldsl Il998t 
[Reynolds fc Keohanelll99^ . 

As was discussed previously, differing lines of sight can 
sample different paths through the ISM, so there is a 
chance for real variations in measurements across 
an SNR. There can also be small systematic variations 
introduced by fitting the regions with the models dis¬ 
cussed above. Each region may be composed of multiple 
plasma components, so fitting to a single plasma model 
may introduce a systematic error or a bias. To show 
this, we performed fits to simulated spectra produced 
with a range of Nh- We produced the simulated spec¬ 
tra using the sedov model in xspec (version 12.8.1 using 
NEIVERS 1.1), which models the total emission from an 
SNR undergoing adiabatic expansion (|Borkowski et al.l 
l2001h . The complete set of simulated spectra cover a 
grid of inputs for shock temperature (Tg = 0.3 —0.9 keV), 
electron temperature (Tes = 0.15 — 0.45 keV), ionization 
age (tq = 1 — 10 X 10^^ s cm“^), and Nn (0.5 — 1.0 x 10^^ 
cm“^). For each parameter set, we simulated 200 spec¬ 
tra and then binned these to have at least 50 counts per 
bin. We fit each spectrum with a pshock model, using 
the tbabs absorption model with the wilms abundance 
table. We present the average of the best fit values for 
each parameter set in Table [TJ along with the standard 
deviation of all fit values. 

The results in Table [T] show that while a single plasma 
model approximately recovers the assumed Vh value 
from a complex plasma, it can introduce significant sta¬ 
tistical errors and suggests a possible bias. In particular, 
the Vh values we measure from the simulated data are 
close to but often lower than the assumed values. This is 
in fact not surprising: fitting a single-temperature model 
to a simulation that is meant to represent the entire 
remnant (i.e., a Sedov spectrum which integrates over 
a range of temperatures) cannot adequately describe the 
spectrum or yield a very accurate value of the hydrogen 
column density. 

To minimize these effects in the analysis of the ac¬ 
tual data, we take a three-pronged approach. First, 
we extract spectra from as small regions of the rem¬ 
nant as possible, to avoid creating complicated, multi¬ 
temperature regions with multiple plasma components. 
Second, we perform fits with thermal, non-thermal, and 
mixed thermal/non-thermal models, to capture the spec¬ 
tral characteristics of the regions correctly. Third, we 
use Bayesian statistical tools (discussed in Section 14.2|) 
to combine measurements from different spatial regions 
as well as from different spectral model results to assess 
any systematic uncertainties in the iVu measurement for 
a given remnant. It is important to note, however, that 
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the simulated spectra results on their own are not the 
sole motivation for the Bayesian analysis; if we believed 
that the simulated results were immediately compara¬ 
ble to the Chandra data results, then we could use the 
simulated results to establish the magnitude of a bias 
term. In Section 5, we will use the total uncertainties 
we determine for the Nn towards each remnant to more 
accurately constrain the N-a/Ay relationship. 

3. CHANDRA DATA PIPELINE AND PROCESSING 

3.1. Pre-Processing 

For each of the selected SNRs, we downloaded archival 
data from the Chandra archive. For remnants with mul¬ 
tiple pointings and exposures, we chose one observation 
to eliminate any need to co-add multiple exposures, and 
placed preference on newer observations that were long 
enough to provide multiple regions with sufficient counts 
for spectral analysis. We ran the data through two 
pipeline scripts that were designed to automate the pro¬ 
cessing procedures. All scripts utilized Cl AO tools used 
version 4.5. Our first script reapplied the latest calibra¬ 
tion and produced a new levels 1 event file. Because the 
majority of remnants were taken with exposure mode 
FAINT or VFAINT, we cleaned the ACIS background 
using procedures that were appropriate for these modes. 
Finally, we filtered the levels I file for bad grades and ap¬ 
plied the good time intervals to generate the new level=2 
event file. 

3.2. Region Selection and Data Processing 

We used the new level=2 event file to select regions 
for spectral extraction. We selected regions using ds9 
and chose areas that had > 10,000 counts, as well as 
sampling a variety of lines-of-sight across the remnant. 
In some remnants, it was preferable or necessary to select 
regions with non-thermal emission. We processed these 
using procedures identical to the thermal regions, but fit 
them with non-thermal models. In general, regions were 
selected using the morphology of the remnant as a guide, 
such that the regions would contain one variety of plasma 
(the importance of which is highlighted in Section 2). We 
extracted a spectrum from each region and then grouped 
it to a minimum of 25 counts per bin. 

As a typical example of the regions selected for a rem¬ 
nant, we show SNR G109.1 — 1.0 in FigurelU with a num¬ 
ber of the selected regions marked with circles and the re¬ 
gions selected for local background shown with a dashed 
circle and rectangle. In the full analysis of G109.1—1.0, 
we used a total of eleven regions. However, for clarity and 
simplicity, we chose to highlight only six in this figure. 
We show the resultant binned spectra from two of these 
regions in Figure[2l again as two representative examples 
of the types of spectra we encountered in the analyses. 
Note that small contributions from model components 
(such as the power-law component in Region 3 of Fig¬ 
ure [2]) may appear insignificant, but are important for 
producing acceptable fits. 

We fit a variety of models (discussed in Section [2]) to 
each grouped spectrum. We fit each model with default 
ISM abundances, but it was often necessary to allow 
some of the plasma abundances to vary in order to prop¬ 
erly fit regions with ejecta-enriched plasma. The num¬ 
ber and variety of the free abundances were allowed to 


change from model to model for a given region, and only 
the plasma models’ abundances were allowed to vary. For 
non-thermal regions, this process was the same, but sim¬ 
plified by the fact that there were no free abundance 
parameters possible for those models. 

Despite our best efforts to select regions that would be 
composed of only one variety of plasma (and thus could 
be fit using one model), sometimes this was infeasible (for 
example, if the regions had to be large in order to contain 
an adequate number of counts). In such cases we needed 
to add other components to the models. Most commonly, 
this required adding a nonthermal model to the thermal 
model to compensate for regions near nonthermal fila¬ 
ments or a central pulsar wind nebula (PWN). However, 
in some cases we added Gaussian features to fit spectral 
lines that were either not properly fit or non-existent in 
the thermal model, or to account for ejecta mixing into 
a non-thermal region, causing emission lines to appear 
on top of a power-law spectrum. For the most part, the 
data were good enough that we could find regions that 
allowed for a single model to produce an acceptable fit, 
but when included, these additional components did not 
make statistically significant changes in the best-fit value 
of A^h- 

We determined if a spectral fit was good both by vi¬ 
sually checking the final fit, as well as by using the 'X^ jv 
fit statistic. We show in Figure |3] an example of an ac¬ 
ceptable fit compared to a poor fit. We did not place a 
hard upper limit on the fit statistic because some large 

jv values were dominated by a single feature that was 
not captured by the model but did not affect the inferred 
A^h- Nevertheless, nearly all of the fits used in our final 
analysis had jv < 2. 

4. RESULTS 

4.1. Results of The Spectral Analysis 

For each SNR region, we included all models that pro¬ 
duced acceptable fits. As a typical example, we show 
in Table [5] the complete fit results for six of the regions 
of G109.I—1.0. It is evident from this table that vnei 
and vpshock models with thawed abundances were com¬ 
monly the best models for remnants with thermal re¬ 
gions, with a powerlaw component added in when neces¬ 
sary. The raymond and mekal models ranged from fairly 
successful to very poor, which presumably reflects sig¬ 
nificant non-equilibrium ionization contributions in the 
data. For SNRs with completely non-thermal regions, 
the powerlaw model was the default model with con¬ 
sistent success (Two of the remnants, G120.1-I-1.4 and 
G04.5-I-6.8, were better fit by an srcut model). All er¬ 
rors presented on the fit A^h values in Table [ 2 ] are 90% 
confidence ranges. The complete fit results for all rem¬ 
nants used in our analysis have been released as a Zenodo 
dataset (DOI: I0.5281/zenodo.I7I83). 

4.2. Bayesian Analysis to Determine A^h and Its 
Uncertainty 

Each of the fits within a region have formal (and often 
asymmetric) uncertainties, as we showed in Table [2] for 
SNR G109.I —1.0. From these fits, the dispersion in the 
measurements arising from spatial sampling or from dif¬ 
ferent continuum models can be calculated. For a num¬ 
ber of remnants, this contribution is comparable to the 
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Table 2 

Example Fit Values for SNR G109.1—0.0 


Region 

Model 

Thawed Elements 

Fit Nn 
(10^^ cm^) 


1 (NE) 

vnei + powerlaw 

Mg, Si 

_o 07 

146.94/97 


vpshock + powerlaw 

Mg, Si, S 

1 nc+0.08 

l.UO_o 07 

143.81/96 

3 (N) 

nei + powerlaw 

None 

U.yO_o 09 

132.71/111 


vnei + powerlaw 

Mg, Si 

0 QQ+O-IO 

u.»o_o 16 

130.55/109 


pshock + powerlaw 

None 

u.yo_o 04 

131.83/111 


vpshock + powerlaw 

Mg, Si 

u.oy_o 06 

126.88/109 

4 (E) 

vnei + powerlaw 

Mg, Si 

u.yu_o 07 

160.49/124 


vpshock + powerlaw 

Mg, Si 

0.931H? 

153.14/124 

6 (Genter) 

vnei 

Mg, Si, S 

o.selH? 

138.73/110 


vpshock 

Mg, Si, S 

0 71+0-08 
'-'•'-‘--0.05 

130.71/110 

9 (SE Inner) 

vnei 

Mg, Si, S 

-L-^4_o 07 

139.21/117 


vpshock 

Mg, Si, S 

1 Qq + 0.05 
i.zo_o 05 

137.94/117 

11 (SE Outer) 

nei 

None 

n Q7+0.06 

t)-of_0.06 

142.52/111 


vnei 

Mg, Si 

U.50_o,o7 

142.15/109 


pshock 

None 

t)-oo_0.06 

130.27/111 


vpshock 

Mg, Si 

U.5y o,o8 

129.88/109 


> 


3 

c 

3 

o 

U 


S 

o 

Z 




1 2 
Energy (keV) 


1 2 
Energy (keV) 


Figure 3. Two example model fits to a spectrum of SNR G109.1-1.0 extracted from Region 6. The acceptable fit (left panel) is a 
vnei model, whereas the unacceptable fit (right panel) is a vraymond model in this particular example. Each model has Mg, Si, and S 
abundances as free fit parameters. 


formal uncertainties. In others, however, the difference 
between the various measurements of A^h is significantly 
larger than the formal uncertainties, pointing to system¬ 
atic uncertainties originating from spatial sampling, from 
the choice of model, or both. For the purposes of our 
analysis, we do not distinguish between the variance of 
the ISM along different lines of sight and the system¬ 
atic error introduced by the incomplete plasma models 
and we combine these different types of syst ematic er¬ 
ror into a single term (see, e.g.. lSinervoll2003l ). Our aim 
in this section is to determine the most likely A^h value 
for each remnant, as well as a measure of the combined 
formal and systematic uncertainties in this quantity. We 
accomplish this by finding the parameters of the under¬ 
lying A^h distribution for each remnant that is consistent 
with our sample of measurements. We use this distribu¬ 
tion to find the most likely value of the hydrogen column 
density and its uncertainty. 


We start with a parametric form of the A^h distribution 
and use the set of measurements from the individually fit¬ 
ted regions in order to estimate its parameters. We take 
the assumed underlying distribution to be a Gaussian 


P(A^h;ct, A^Hc) = Cexp 


(iVn-iVnc)'] 

2a2 


( 1 ) 


with a mean A^hc and a standard deviation a that can 
be different for each remnant. In this and the following 
expressions, C is a proper normalization constant such 
that 

P(A^h; O', A^Hc)dAfH = 1- (2) 

We also need to model the individual measurement un¬ 
certainties where i represents a particular SNR 

region/model combination that yields a single measure¬ 
ment of A^h for that region. In general, the surface for 
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Table 3 

Measured Column Densities for All SNRs 


SNR 

Number of 
Regions/Models 

Ah 

(10^2 

Error 

cm^) 

GO.0+0.0 

3 

9.86 

0.34 

G04.5-I-6.8 

4 

0.54 

0.011 

G06.4-0.1 

22 

0.63 

0.19 

G53.6-2.2 

2 

0.67 

0.03 

G54.1-I-0.3 

1 

2.55 

0.04 

G69.0-I-2.7 

1 

0.45 

0.02 

G109.1-1.0 

30 

0.90 

0.14 

Gill.7-2.1 

8 

1.17 

0.16 

G116.9+0.2 

2 

0.92 

0.07 

G119.5+10.2 

1 

0.38 

0.11 

G120.1+1.4 

14 

0.80 

0.10 

G130.7+3.1 

6 

0.54 

0.01 

G184.6-5.8 

1 

0.30 

0.02 

G260.4-3.4 

3 

0.35 

0.15 

G263.2-3.3 

1 

0.03 

0.01 

G327.6+14.6 

2 

0.16 

0.01 

G332.4-0.4 

9 

0.87 

0.35 


each measurement is not simple or symmetric around the 
minimum, leading to asymmetric formal errors in these 
individual measurements. However, around the mini¬ 
mum, it is accurate to represent the likelihood using two 
half Gaussians with different dispersions that smoothly 
connect at the most likely column density Nno^i for each 
measurement; i.e., 


Pi(data|fVH) = 



(Ah— iVHo.i)^ 


(Ah—Who, i)^ 

-- 

_L NTH 1 


Nh < fVno.i 
.^H > .^HO.i 


( 3 ) 


where “data” stands for the most likely column density 
Nno.i and the two associated uncertainties and 

for the i-th SNR region/model combination. 

We want to calculate the quantity P(tT, A^hcI data), 
which measures the posterior likelihood of the parame¬ 
ters of the Wh distribution, given the observations. Using 
Bayes’ theorem, we can write this as 


P(cr, A^Hcldata) = C' 2 P(data|cr, 7VHc)P(cr)P(A^Hc) , (4) 


where C 2 is a normalization constant and P(cr) and 
-P(-^Hc) are the priors over the values of the Gaus¬ 
sian dispersion, cr, and the peak of the IVh distribu¬ 
tion, IVhc- Here, “data” stands for the ensemble of the 
NHo,i, and a+^NH,i values for a particular SNR. 

We take a flat prior over the Gaussian dispersion a be¬ 
tween (Tmin that is equal to 0.1 times the smallest formal 
uncertainty obtained from a spectral fit for each remnant 
and Umax that is equal to 10 times the largest difference 
between two IVh measurements for each remnant: 


P(a) 


0 , 

1 

*^max <^min ^ 


0 , 


O' < CTmin 

<^min cr <C CTjjjax 

O > ■ 


( 5 ) 


Similarly, we take a flat prior over the centroid of the iVn 
distribution TVhc that spans the range from 0.1 times the 
smallest A^h measurement to 10 times the largest TVh 
measurement per source. These limits ensure that the 
particular minimum and maximum values of the prior 
distribution do not affect the results. 


In equation (jH), the quantity P (data | cr, A^hc) measures 
the likelihood that we will make a particular set of mea¬ 
surements for the column density given the values of the 
parameters of the column density distribution. We need 
to estimate this quantity, given the measurement likeli¬ 
hoods given in equation ([3]). We will assume that each 
measurement is independent, so that 

P(data|cr, A^hc) = 

f dAtHPi(data|A^//)P(A^H; O', A^hc) • (6) 

i 

Combining this last equation with equation (jd]), we ob¬ 
tain the posterior likelihood 

P(cr, AtRcldata) = C'P(cr)P(A^Hc)x 

Y[ j dA^HP,(data|A^H)P(A^H;^,AtHc) , (7) 

i 

where C is the overall normalization constant. We can 
use equation [7] to determine the parameters a and A^hc 
given the individual region/model fits for each remnant. 
The dispersion (tr) is then a measure of the systematic 
uncertainty associated with the remnant. If a remnant 
has a significant systematic uncertainty, that contribu¬ 
tion is considered when we define the most likely value 
of the A^h for the SNR and its uncertainty. 

In Figure |H we show as examples the Bayesian cred¬ 
ible regions over the parameters of the A^h distribution 
that correspond to the measurements for G109.1 — 1.0 and 
G4.5+6.8. We show with a dot the most likely values of 
the centroid and dispersion; i.e, the peak of the posterior 
likelihood given by equation ([7]). In the case of G4.5-1-6.8, 
the best-fit value of the dispersion is very close to zero 
and this parameter is consistent with being zero within 
the p = 0.05 credible region, indicating a negligible level 
of systematic uncertainty. On the other hand, for rem¬ 
nants such as G109.1—1.0, the best fit value of the dis¬ 
persion is not consistent with being zero, indicating the 
presence of systematic uncertainties arising from spatial 
sampling across the remnant or the differences among 
spectral models. 

The previous Bayesian analysis to quantify the system¬ 
atic error term yields two categories of remnants. In the 
first category the systematic dispersion {a) is consistent 
with zero, so we assign a zero systematic uncertainty, 
i.e., set P{N}i;a, Nuc) = HNn — A^hc) in equation (|6l). 
This allows us to find a properly weighted average of the 
A^h measurements using their formal uncertainties and 
report the best-fit value and its formal uncertainty in Ta¬ 
ble m In the second category of remnants, we compute 
the posterior likelihood over their A^h by weighing each 
Gaussian distribution with a given centroid Nhc and dis¬ 
persion (T with the likelihood calculated in equation 0 
that those pair of parameters represent the observed A^h 
values: 

P{Nii) = 

J J P(AfH;cr, AfHc)H((7, A^Hc|data)dAfHcdcr (8) 

We use this distribution to infer the best-fit values of 
A^h and its uncertainty for each source and report this in 
Table [3 
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Figure 4. p = 0.32 and p = 0.05 Bayesian credible regions for the central value and the dispersion of an underlying Gaussian distribution of 
A^h measurements for SNRs G109.1-1.0 (left) and G04.5-I-6.8 (right). The dispersion of the distribution reflects the systematic uncertainty 
in the measurements. The dispersion contours for G04.5-I-6.8 are consistent with 0.0, which indicates that systematic uncertainties are not 
significant. For G109.1—1.0, the most likely dispersion is not consistent with zero, indicating a level of systematic error comparable to the 
formal errors. These are included in the uncertainties in Table [3| following the discussion in Section 4. 


5. APPLICATION TO THE Nu/Ay RELATION 

We presented in the previous section the hydrogen col¬ 
umn densities toward a large sample of SNRs that we 
measured using the Chandra data archive. We now com¬ 
bine these Vjj measurements with the optical extinction 
measurements presented in iGuver fc Ozell (|2009f ) for our 
sample of SNRs to determine a relation between these 
two quantities that accounts for systematic errors. The 
optical extinction toward the remnants in our sample was 
determined through several different methods, which we 
list in Table m Most of these methods involve measuring 
the intensity ratio of two emission lines, for which the in¬ 
trinsic ratio is known. Comparing the observed intensity 
ratio to the intrinsic ratio yields a reddening, which is 
then converted into a measurement of the optical extinc¬ 
tion. These methods are challenging, in general, because 
a high signal-to-noise ratio spectrum is necessary to ob¬ 
tain the intensity ratios. 

Several emission line pairs are frequently used for this 
purpose. Ha (6563 A) to Hp (4861 A) line ratio, re¬ 
ferred to as the Balmer decrement method, is one of the 
most well known and reliable ones among these pairs. 
Othe r emission line ratios from th e SII multiplet ([Milled 
Il968f) and FE[II] IR transitions (jOliva et al.l 1198^ are 
also used in a similar way. For all of these pairs of lines, 
the ratio depends only very weakly on the temperature 
and density of the emitting plasma, leading to minimal 
uncertainties in the calculation of the ratio (jOsterbrockl 
[TM^ILeaueil^lMi^ . 

Another common method is to use stars near an SNR 
to estimate the extinction toward that remnant. This 
method relies on having an existing estimate of the dis¬ 
tance to the remnant, and then identifying stars at sim¬ 
ilar distances. If suitable stars can be observed, the ex¬ 
tinction meas urements of those stars can be appli ed to 
the remnant (|Koo et al.ll2008t rR,uiz-LaDuentell20r)4D . 

In Figure [HI we plot the hydrogen column density mea¬ 
surements and their uncertainties presented in Table |3] 
against the measurements of the optical extinction sum¬ 


marized in Tabled! Following IGiiver fc Ozell l|2009f) . we 
assign a 15% error to the optical extinction measure¬ 
ments for which uncertainties have not been reported 
(denoted by in Tabled]). In addition, upon inspecting 
the uncertainties of individual measurements in Table dl 
we note that a best-fit line between these two quantities 
will be heavily influenced by a small number of remnants 
where Vh measurements have very small (< 5%) formal 
uncertainties. Even though the sample of regions and 
the range of models we considered in the spectral fits 


T able 4 _ 

Ay Values from IGiiver Ik Ozell II2009I ) 


SNR 

Av 

(mag) 

Error“ 

(mag) 

Method 

Ref. 

GO.0+0.0 

29 

2 

Nearby Stars 

(1) 

G04.5-I-6.8 

2.5 

0.9 

Fell Ratio 

(2) 

G06.4-0.1 

3.57 

0.47 

H„/H^ 

(3) 

G53.6-2.2 

3.57 

0.47 

Hc/H^ 

(3) 

G54.1-I-0.3 

8.0 

0.70 

Nearby Stars 

(4) 

G69.0-I-2.7 

2.48 

- 

Hc/H^ 

(5) 

G109.1-1.0 

3.15 

0.65 

Hc/H^ 

(6) 

Gill.7-2.1 

5.0 

0.40 

SII ratio 

(7) 

G116.9+0.2 

2.70 

0.5 

Hc/H^ 

(8) 

G119.5+10.2 

1.27 

0.41 

Extinction Map 

(9) 

G120.1+1.4 

1.86 

0.12 

Nearby Stars 

(10) 

G130.7+3.1 

2.11 

- 

Nearby Stars 

(11) 

G184.6-5.8 

1.55 

0.186 

Lya Absorption 

(12) 

G260.4-3.4 

2.60 

- 

Nearby Stars 

(13) 

G263.2-3.3 

0.38 

- 

Hc/H^ 

(14) 

G327.6+14.6 

0.34 

- 

HI/GC 

(15) 

G332.4-0.4 

4.70 

0.90 

Fell Ratio 

(2) 


“Errors desiRn ated as ’ are taken a s 15%. _ 

References: fU IRredehl fc 'fruemperl 1119941 k (2) lUliva et al.l 
179891): ('airLong et al.l (ITMII') : GilKoo et al.l OOOSi'): 1 5') 
HesteiyfclA]llgjni[J []A89l ) : (6ilFesen fc^urfor d 119951 ): (7) 
Hurford fc FesenI ll 1999): (8 1 IFesen et al.l Ill9 9'm: (9) 
Mavromatalas_et^ I ll200Qlb (10) IRuiz-Lapu*^^ feooll: (11) 

al.l II1988I): ISoller man et aTT 20001): fl3) [GorensteinI 

(119751) : (14) IManchester et al.llll97^) : flS HRaymond et al.l (119951) . 
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Figure 5. The A^h and Ay measurements for the SNRs in the Chandra archive along with the best fit linear model, which gives 
A^h = 2.87±0.12 X 10^^ Ay. The A^h measurements are presented in Table[3]and have a minimum error of 5%, while the Ay measurements 
(listed in Table l4t are taken from IGuver fc Ozell II2009I V 


did not allow us to measure a systematic uncertainty for 
these remnants, considering the many possible sources of 
systematic uncertainty that are usually present in the de¬ 
termination of the hydrogen column densities, we assign 
a 5% error to the iVn measurement of SNR G130.7-1-3.1, 
SNR G69.0 -b 2.7, SNR G004.5 -b 06.8, SNR G53.6 - 2.2, 
SNR G54.1 -b 0.3, and SNRG 00.0 -b 0.0 to avoid biasing 
the results. 

Using these data, we obtain a best-fit linear relation 
between the hydrogen column density and the optical 
extinction that is described by 

(Vh = 2.87 ± 0.12 X 10^1 Ay cm-2 (9) 

where Ay is in magnitudes and the error represents the 
1-cr uncertainty. We show the best-fit line in Figure [S] 

We can test the effect of imposing a floor to the Nh 
uncertainties by also fitting to the raw values. This in¬ 
creases the best fit value of the N-a/Ay to (2.92±0.11) x 
10^^ cm“^, which is within 1-cr of the best fit value pre¬ 
sented above, but dominated by remnants with small er¬ 


rors. Finally, we note that the remnants with either very 
low or high Nn will have the greatest leverage on the 
final fit. To test the magnitude of this effect, we remove 
the points with the highest and lowest values of the hy¬ 
drogen column density. This results in a best fit value of 
(2.76±0.13) X 10^^ cm“^, again within 1-a of the best fit 
value using all of the data and a minimum error imposed 
on the A^h measurements. 

These results indicate that a consistently analyzed 
sample of SNRs using ISM abundances for spectral fit¬ 
ting produces an N-^/Ay relation that is significantly 
higher than previous estimates. This has the effect of 
increasing the iVu derived from existing E(B-V) mea¬ 
surements. Conversely, using the new relation to esti¬ 
mate the brightness of an optical counterpart of an X- 
ray source will result in less extinction compared to es¬ 
timates calculated with previous values of the Nn/Ay 
relation. As antici pated by Giiver fc Ozell (I2009D and 
further discussed by IWatsonI ( 20111) . the primary reason 
for this difference is the change in the abundances that 
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are used for the ISM. The default metallicity library of 
the frequently used spectral analysis package xsp e c uses 
the solar abundances given by I Anders fc Grevessd (|1989(1 
and this is indeed the main library that is utilized in the 
studies of the spectral properties of X-ray sources. How¬ 
ever, these values are known to be on average 45% higher 


a,g e 4070 1 

iTdM) 


than the updated values bv lAsplund et al.l II^UUt)l) using 
the solar spectrum and bv iWilms et al.l ( 200011 for the 
ISM and, therefore, result in a smaller A^h value for a 
given amount of total absorption due to the matter in 
the ISM. We showed here in our systematic analysis that 
a change in the abundance table results in Ri 30% change 
in the coefficient of the linear relation between the hy¬ 
drogen column density and optical extinction. 

As a final caveat, we note that when the hydrogen col¬ 
umn density is measured using the ISM abundances (e.g., 
from IWilms et al.l l2000f) , the relation we presented here 
should be used to predict or compare with the extinction 
in the optical band. On the other ha nd, if an TVh is found 
using solar abundances (e.g., from lAnders fc Grevessi 
Il989f) , then the relation reported in lGiiver fc Ozell (12009( 1 
should be employed for self consistency. 


populations of sources that are bright in both optical and 
X-ray bands, such as blazars. The measurement of the 
hydrogen column density and the optical extinction in 
such a population would provide an independent mea¬ 
surement of the relation between these two quantities 
and help cross check the results. 
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